Rによる回帰分析 (Regression using R)

回帰 (Regression)

以下のYouTube動画を視聴して、回帰 (Regression)の概念を理解してください。動画は英語ですが、日本語字幕機能をオンにすることもできます。 Please watch the following youtube video to get the idea of regression analysis

• Selection Bias: Will You Make More Going to a Private University?

実証例のデ一タ可視化 (Data Visualization)

研究課題：少人数教育が学生の成績に与える定量的効果は何か？

以下のデータセットCASchoolsを使用します。 ここで使用するデータは、カリフォルニア州の420のK-6とK-8の全地区のデータで、1998年と1999年のデータが利用可能です。テストのスコアは、5年生を対象に実施されたスタンフォード9標準化テストのものです。 学校の特徴（地区全体の平均値）には、入学者数、教師数（「フルタイム換算」で測定）、教室あたりのコンピュータの数、生徒一人あたりの支出が含まれます。 生徒の人口統計学的変数は、地区全体で平均化されています。 人口に関連する変数には、生活保護プログラムであるCalWorks（旧AFDC）の生徒の割合、割引価格の昼食の対象となる生徒の割合、英語学習者（英語が第二言語である生徒）の割合が含まれています。

We are going to use the following dataset CASchools. The data used here are from all 420 K-6 and K-8 districts in California with data available for 1998 and 1999. Test scores are on the Stanford 9 standardized test administered to 5th grade students. School characteristics (averaged across the district) include enrollment, number of teachers (measured as “full-time equivalents”, number of computers per classroom, and expenditures per student. Demographic variables for the students are averaged across the district. The demographic variables include the percentage of students in the public assistance program CalWorks (formerly AFDC), the percentage of students that qualify for a reduced price lunch, and the percentage of students that are English learners (that is, students for whom English is a second language).

# データはAERというパッケージに含まれています
# 下記のパッケ一ジをlibrary()関数で読み出します
library(AER)

# 利用したいパッケ一ジがなかったら、
# 下記のようにパッケ一ジをインスト一ルしておきます (初回のみ); 例として:
# install.packages("AER")

# デ一タを読み込む
# AER::CASchools 
# load the data set onto the workspace
data("CASchools") 

# dataframe structure
str(CASchools)

'data.frame':   420 obs. of  14 variables:
 $ district   : chr  "75119" "61499" "61549" "61457" ...
 $ school     : chr  "Sunol Glen Unified" "Manzanita Elementary" "Thermalito Union Elementary" "Golden Feather Union Elementary" ...
 $ county     : Factor w/ 45 levels "Alameda","Butte",..: 1 2 2 2 2 6 29 11 6 25 ...
 $ grades     : Factor w/ 2 levels "KK-06","KK-08": 2 2 2 2 2 2 2 2 2 1 ...
 $ students   : num  195 240 1550 243 1335 ...
 $ teachers   : num  10.9 11.1 82.9 14 71.5 ...
 $ calworks   : num  0.51 15.42 55.03 36.48 33.11 ...
 $ lunch      : num  2.04 47.92 76.32 77.05 78.43 ...
 $ computer   : num  67 101 169 85 171 25 28 66 35 0 ...
 $ expenditure: num  6385 5099 5502 7102 5236 ...
 $ income     : num  22.69 9.82 8.98 8.98 9.08 ...
 $ english    : num  0 4.58 30 0 13.86 ...
 $ read       : num  692 660 636 652 642 ...
 $ math       : num  690 662 651 644 640 ...

# 下記いくつの データサイエンス パッケ一ジを library()関数で読み出します
# data science packages
library(tidyverse)
library(dplyr)
library(magrittr)

# 少人数教育を測る変数(st.ratio: student-teacher ratio)を作る
# compute STR and append it to CASchools
# dplyr::mutate
CASchools %<>% mutate(st.ratio = students/teachers)

# クラスサイズ変数の構築: d
# construct the class-size variable: d
# d =small_class, STR <20; d =large_class, STR>=20
CASchools %<>% mutate(d =rep("large_class", length(CASchools$district)))
CASchools$d[CASchools$st.ratio < 20] <- "small_class"

# ｄをカテゴリ変数にする
# make variable d into a categorical variable
CASchools$d %<>% as.factor

# 学生のパフォーマンスを測定するアウトカム変数 `score` を作成する
# score = 読解力と数学のテストの平均点です
# compute TestScore and append it to CASchools
CASchools %<>% mutate(score = (read+math)/2)

# dataframe structure
str(CASchools)

'data.frame':   420 obs. of  17 variables:
 $ district   : chr  "75119" "61499" "61549" "61457" ...
 $ school     : chr  "Sunol Glen Unified" "Manzanita Elementary" "Thermalito Union Elementary" "Golden Feather Union Elementary" ...
 $ county     : Factor w/ 45 levels "Alameda","Butte",..: 1 2 2 2 2 6 29 11 6 25 ...
 $ grades     : Factor w/ 2 levels "KK-06","KK-08": 2 2 2 2 2 2 2 2 2 1 ...
 $ students   : num  195 240 1550 243 1335 ...
 $ teachers   : num  10.9 11.1 82.9 14 71.5 ...
 $ calworks   : num  0.51 15.42 55.03 36.48 33.11 ...
 $ lunch      : num  2.04 47.92 76.32 77.05 78.43 ...
 $ computer   : num  67 101 169 85 171 25 28 66 35 0 ...
 $ expenditure: num  6385 5099 5502 7102 5236 ...
 $ income     : num  22.69 9.82 8.98 8.98 9.08 ...
 $ english    : num  0 4.58 30 0 13.86 ...
 $ read       : num  692 660 636 652 642 ...
 $ math       : num  690 662 651 644 640 ...
 $ st.ratio   : num  17.9 21.5 18.7 17.4 18.7 ...
 $ d          : Factor w/ 2 levels "large_class",..: 2 1 2 2 2 1 2 1 2 1 ...
 $ score      : num  691 661 644 648 641 ...

# summary statistics
summary(CASchools)

   district            school                  county      grades   
 Length:420         Length:420         Sonoma     : 29   KK-06: 61  
 Class :character   Class :character   Kern       : 27   KK-08:359  
 Mode  :character   Mode  :character   Los Angeles: 27              
                                       Tulare     : 24              
                                       San Diego  : 21              
                                       Santa Clara: 20              
                                       (Other)    :272              
    students          teachers          calworks          lunch       
 Min.   :   81.0   Min.   :   4.85   Min.   : 0.000   Min.   :  0.00  
 1st Qu.:  379.0   1st Qu.:  19.66   1st Qu.: 4.395   1st Qu.: 23.28  
 Median :  950.5   Median :  48.56   Median :10.520   Median : 41.75  
 Mean   : 2628.8   Mean   : 129.07   Mean   :13.246   Mean   : 44.71  
 3rd Qu.: 3008.0   3rd Qu.: 146.35   3rd Qu.:18.981   3rd Qu.: 66.86  
 Max.   :27176.0   Max.   :1429.00   Max.   :78.994   Max.   :100.00  
                                                                      
    computer       expenditure       income          english      
 Min.   :   0.0   Min.   :3926   Min.   : 5.335   Min.   : 0.000  
 1st Qu.:  46.0   1st Qu.:4906   1st Qu.:10.639   1st Qu.: 1.941  
 Median : 117.5   Median :5215   Median :13.728   Median : 8.778  
 Mean   : 303.4   Mean   :5312   Mean   :15.317   Mean   :15.768  
 3rd Qu.: 375.2   3rd Qu.:5601   3rd Qu.:17.629   3rd Qu.:22.970  
 Max.   :3324.0   Max.   :7712   Max.   :55.328   Max.   :85.540  
                                                                  
      read            math          st.ratio               d      
 Min.   :604.5   Min.   :605.4   Min.   :14.00   large_class:181  
 1st Qu.:640.4   1st Qu.:639.4   1st Qu.:18.58   small_class:239  
 Median :655.8   Median :652.5   Median :19.72                    
 Mean   :655.0   Mean   :653.3   Mean   :19.64                    
 3rd Qu.:668.7   3rd Qu.:665.9   3rd Qu.:20.87                    
 Max.   :704.0   Max.   :709.5   Max.   :25.80                    
                                                                  
     score      
 Min.   :605.5  
 1st Qu.:640.0  
 Median :654.5  
 Mean   :654.2  
 3rd Qu.:666.7  
 Max.   :706.8  
                

# 数値変数を引き出し、その相関係数を計算します。
# これは、どのような変数を回帰式に含めるべきかのヒントを与えてくれる。
# 次のようなデータ可視化により、これらの変数間の関係を簡単に識別することができる。

# (Pearson) linear correlation: outcome, target, controls
target_controls <- select(CASchools, score, st.ratio, calworks:english) 
cor_matrix <- cor(target_controls)

# Kendall coefficient: dependence
# cor(target_controls, method = c("kendall"))

# correlation matrix plot
# install.packages("corrplot")
library(corrplot)
corrplot::corrplot(cor_matrix, method = "circle")

#corrplot::corrplot(cor_matrix, method = "color", addCoef.col = "grey")

研究課題：少人数教育 (st.ratio) が学生の成績 (score) に与える定量的効果は何か？

# キレイなグラフを作るためのパッケージ
# data visualization, # qqplot2::qqplot
library(ggplot2)

# scatter plot
#install.packages(c("ggExtra", "ggthemes"))
library(ggExtra)
library(ggthemes)

p <- ggplot(data =CASchools, aes(y =score, x =st.ratio)) +
  geom_point(alpha = .5) +
  theme_tufte(ticks=F)
ggMarginal(p, type = "histogram", fill="transparent")

データ可視化を上手に行うためのお役立ちリンク集 (Useful links for doing good data visualization)

• Schwabish J.A. (2014), “An Economist’s Guide to Visualizing Data,” Journal of Economic Perspectives, Vol. 28, 209–234.

• Data Visualization – a practical introduction

• Tufte in R – The visual display of quantitative information

トップへ戻る

Rによる回帰分析 (Regression in R)

単回帰分析 (Simple regression illustrated through interactive graphs)

• リンク: Estimating the coefficients of a simple linear regression (2D plot)

• リンク: Adding estimated hyperplane in the 3D plot

 

回帰式

\[ \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i \]

lm(y ~ x)

回帰分析は lm() を利用することで実行できます。 Rで回帰分析を行うと、推定されたバラメ一タ \((\hat{\beta}_0, \ \hat{\beta}_1)\) の標準誤差とともに、有意差検定 (t検定) の結果を自動で出力してくれます。 回帰分析の結果を保存したオブジェクトを summary() で要約すると、標準誤差は standard.error としてレポ一トされており、有意差検定の結果は t-statistics、p-value としてそれぞれの変数に対して表示されています。

処置変数 d (カテゴリ変数)

\(d=1\) if st.ratio<20; otherwise, \(d=0\).

回帰式

\[ \hat{score}_i = \hat{\beta}_0 + \hat{\beta}_1 d_i \]

# regress score on (d: discrete r.v.)
CASchools %$% lm(score ~ (st.ratio<20)) %>% summary

Call:
lm(formula = score ~ (st.ratio < 20))

Residuals:
    Min      1Q  Median      3Q     Max 
-50.496 -14.029  -0.346  12.884  49.504 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        650.077      1.393 466.666  < 2e-16 ***
st.ratio < 20TRUE    7.169      1.847   3.882  0.00012 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 18.74 on 418 degrees of freedom
Multiple R-squared:  0.0348,    Adjusted R-squared:  0.0325 
F-statistic: 15.07 on 1 and 418 DF,  p-value: 0.0001202

クイズ: \(\hat{\beta}_0\) と \(\hat{\beta}_1\) の解釈は？ゼロとは有意に違うのでしょうか？

 

# load the color palettes 
library(wesanderson)
# names(wes_palettes)
# e.g., Moonrise2, GrandBudapest1, GrandBudapest2

# グラフ
# small class size (st.ratio < 20) v.s. large calss size 
# scatter plots, linear specification
sp <- CASchools %>%
ggplot(aes(y =score, x =st.ratio, color =d)) +
  geom_point(alpha = .7) +
  stat_smooth(method = "lm", formula = y ~ x)

# d is a facotr variable with 2 levels; therefore n=2 
sp + scale_color_manual(values=wes_palette(n=2, name="Moonrise2"))

処置変数 st.ratio (連続変数)

回帰式

\[ \hat{score}_i = \hat{\beta}_0 + \hat{\beta}_1 st.ratio_i \]

# regress score on (st.ratio: continuous r.v.)
lm.fit.1 <- CASchools %$% lm(score ~ st.ratio)
summary(lm.fit.1)

Call:
lm(formula = score ~ st.ratio)

Residuals:
    Min      1Q  Median      3Q     Max 
-47.727 -14.251   0.483  12.822  48.540 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 698.9329     9.4675  73.825  < 2e-16 ***
st.ratio     -2.2798     0.4798  -4.751 2.78e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 18.58 on 418 degrees of freedom
Multiple R-squared:  0.05124,   Adjusted R-squared:  0.04897 
F-statistic: 22.58 on 1 and 418 DF,  p-value: 2.783e-06

クイズ: \(\hat{\beta}_0\) と \(\hat{\beta}_1\) の解釈は？ゼロとは有意に違うのでしょうか？

 

# グラフ
# 網掛け部分は95%信頼区間を表します
# scatter plot, linear specification
CASchools %>%
ggplot(aes(y =score, x =st.ratio)) +
  geom_point(alpha = .5) +
  stat_smooth(method = "lm", formula = y ~ x, color='#D55E00')

トップへ戻る

重回帰分析 (Multiple regression)

次の2つのグラフは、なぜ共変量 (control variable) english を含めるべきかを示しています。 したがって、重回帰分析を利用することで、より選択バイアスの影響が少ない分析ができます。

# motivation: multiple regressions, 2D scatter plot

# english: the share of English learning students
CASchools %>%
ggplot(aes(y =score, x =st.ratio, color =english)) +
  geom_point(alpha = .5) 

# motivation: multiple regressions, 3D scatter plot
# english: the share of English learning students

#install.packages("plotly")
library(plotly)

plot_interactive <- plot_ly(CASchools, x = ~st.ratio, y = ~english, z = ~score,
             marker = list(color = ~english, colorscale = c('#FFE1A1', '#683531'), showscale = TRUE)) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'st.ratio'),
                      yaxis = list(title = 'PctEnglish'),
                      zaxis = list(title = 'score')),
         annotations = list(
           x = 1.13,
           y = 1.05,
           text = 'share of English learning students',
           xref = 'paper',
           yref = 'paper',
           showarrow = FALSE
         ))  

plot_interactive

回帰式

\[ \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_{1i} + \hat{\beta}_2 x_{2i} \]

# 重回帰分析
# multiple regression
lm.fit.2 <- CASchools %$% lm(score ~ st.ratio + english)
summary(lm.fit.2)

Call:
lm(formula = score ~ st.ratio + english)

Residuals:
    Min      1Q  Median      3Q     Max 
-48.845 -10.240  -0.308   9.815  43.461 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 686.03224    7.41131  92.566  < 2e-16 ***
st.ratio     -1.10130    0.38028  -2.896  0.00398 ** 
english      -0.64978    0.03934 -16.516  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.46 on 417 degrees of freedom
Multiple R-squared:  0.4264,    Adjusted R-squared:  0.4237 
F-statistic:   155 on 2 and 417 DF,  p-value: < 2.2e-16

クイズ: なぜ、st.ratio の推定効果が -2.2 (単回帰) から -1.1 (重回帰) に減少します?

不均一分散ロバスト標準誤差 (Heteroskedasticity-robust standard errors)

不均一分散とは?

テストの得点の分散が、生徒と教師の比率とともに増加することが観察されます。

 

estimatrは、一般的によく使われる線形推定量を提供するRパッケージであり、速度と使いやすさを考慮して設計されています。ユーザは、クラスタランダム化、ブロックランダム化、ブロックおよびクラスタランダム化を反映した線形推定量を選択することができる。

estimatr is an R package providing a range of commonly-used linear estimators, designed for speed and for ease-of-use. Users can choose an estimator to reflect cluster-randomized, block-randomized, and block-and-cluster-randomized designs. https://github.com/DeclareDesign/estimatr

se_type: (ロバスト標準誤差の設定)

• classical: homoskedasticity (均一分散標準誤差, lm() 関数で使用される既定の設定です)
• HC0: \(\omega_i = \hat{u}_i^2\). White (1980)
• stata: same as vec(robust) used by Stata
• HC1: \(\omega_i = \frac{n}{n-k}\hat{u}_i^2\)
• HC2: \(\omega_i=\frac{1}{1-h_i}\hat{u}_i^2\)
• HC3: \(\omega_i=\frac{1}{(1-h_i)^2}\hat{u}_i^2\)

 

lm_robust(y ~ x, se_type = "HC2")

library(estimatr)
# estimatr::lm_robust

lm.fit.2.robust <- CASchools %$% 
  lm_robust(score ~ st.ratio + english, se_type = "HC2")
summary(lm.fit.2.robust)

Call:
lm_robust(formula = score ~ st.ratio + english, se_type = "HC2")

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value   Pr(>|t|) CI Lower CI Upper  DF
(Intercept) 686.0322    8.75433  78.365 1.278e-251 668.8241 703.2404 417
st.ratio     -1.1013    0.43417  -2.537  1.156e-02  -1.9547  -0.2479 417
english      -0.6498    0.03111 -20.888  7.919e-67  -0.7109  -0.5886 417

Multiple R-squared:  0.4264 ,   Adjusted R-squared:  0.4237 
F-statistic: 222.8 on 2 and 417 DF,  p-value: < 2.2e-16

st.ratioに対応する不均一分散ロバスト標準誤差は 0.43417 であり、これは均一分散の下で計算された標準誤差 (0.38028) より大きいです。

トップへ戻る

ノンパラメトリック・ブートストラップ (Nonparametric bootstrap)

ブートストラップは、サンプリング分布を構築するための計算アルゴリズムです。The bootstrap is a computational algorithm for constructing sampling distributions.

• 理論的な近似に依存するのではなく、現在のサンプル（経験的データ分布）から再サンプリングを行うことで、サンプリング分布を模倣します。Instead of relying on theoretical approximation, the bootstrap uses resampling from your current sample (empirical data distribution) to mimic the sampling distribution.

 

---

ALGORITHM 1.  The Nonparametric Boostrap

Given data \(\{z_i \}_{i=1}^n := \{y_i, x_i \}_{i=1}^n\), for \(b =1, \ldots, B\):

• Resample with-replacement \(n\) observations \(\{z_i^b \}_{i=1}^n\)

• Calculate your estimate, say, \(\hat{\beta}_b\), using this resampled set.

Then \(\left\{\hat{\beta}_b \right\}_{b=1}^B\) is an approximation to the sampling distribution for \(\hat{\beta}\).

• The standard deviation of your sampling distribution can be approximated by the standard deviation of this bootstrap sample: \[ \text{se}(\hat{\beta}) \approx \text{sd}(\hat{\beta}_b) = \sqrt{\frac{1}{B}\sum_b\left(\hat{\beta}_b -\hat{\beta}\right)^2}. \]

• So long as the estimator is unbiased, which means that \(\mathbb{E}[\hat{\beta}]=\beta\), then we can use this standard error to build the usual \(95\%\) confidence interval: \[ \beta\in \hat{\beta} \pm 2\text{sd}(\hat{\beta}_b). \]

---

 

The folloiwng algorithm applies to the case where \(\mathbb{E}[\hat{\beta}]\neq\beta\).

ALGORITHM 2.  Nonparametric Bootstrap for Confidence Intervals

Given data \(\{z_i \}_{i=1}^n\) and full sample estimate \(\hat{\beta}\) for parameter \(\beta\), for \(b=1, \ldots, B\):

• Resample with-replacement \(n\) observation \(\{z_i^b \}_{i=1}^n\).

• Calculate your estimate, say, \(\hat{\beta}_b\), using this resampled set.

• Calculate the error, say, \(e_b = \hat{\beta}_b - \hat{\beta}\).

Then \(\{e_b\}_{b=1}^B\) is an approximation to the distribution of errors between estimates and their targe.

• To get the \(95\%\) confidence interval for the true \(\beta\), we calculate the \(2.5\)th and \(97.5\)th percentiles of \(\{e_b\}_{b=1}^B\) – say, \(t_{0.025}\) and \(t_{0.975}\) – and set our interval as \[ \left[\hat{\beta} -t_{0.975}, \ \hat{\beta} -t_{0.025} \right]. \]

 

ブートストラップが破綻するのはいつか？ When does the bootstrap break?

• ブートストラップが破綻するのは、経験的データ分布が母集団の分布を適切に近似できていない場合や、対象とする統計量が母集団全体の分布に関する多くの情報を必要とする場合です。If the empirical data distribution is a poor approximation to the population or if the statistic you are targeting requires a lot of information about the entire distribution.

• 例えば、高次元のパラメーターの場合、母集団内の全ての変数間の依存関係を含む結合サンプリング分布を考慮する必要があるため、ブートストラップを信用しない傾向があります。e.g., we tend to not trust bootstrap for high-dimensional parameters because the joint sampling distribution is a function of all the population cross-variable dependencies.

• このような場合には、ブートストラップの代替手法を検討する必要があります。When this happens, we look to alternative versions of the bootstrap.

Bootstrap in action, Take 1

B <- 1000
betas <- c()
for (b in 1:B) {
  samp_b <- sample.int(nrow(CASchools), replace=TRUE)
  reg_b <- lm(score ~ st.ratio + english, data=CASchools[samp_b,])
  betas <- rbind(betas, coef(reg_b))
}

# note that betas is an array
# display the first five \hat{\beta}_b
betas[1:5,]

     (Intercept)  st.ratio    english
[1,]    685.2681 -1.083098 -0.6348671
[2,]    694.6744 -1.476707 -0.6560710
[3,]    687.2626 -1.133849 -0.6610360
[4,]    686.2406 -1.166699 -0.6136197
[5,]    685.7927 -1.108954 -0.6408945

# calculate the standard deviation of this boostrap sample
# apply(matrix, dimcode, function)
apply(betas[, 1:3], 2, sd)

(Intercept)    st.ratio     english 
 8.88565009  0.44120502  0.03180366 

Bootstrap in action, Take 2

# Nonparametric bootstrapping
library(boot)

# create a function producing beta_hat based on resampled observations
# this user-defined function will be used by boot::boot

beta_hat <- function(formula, data, indices) {
  d.resampled <- data[indices,] # allows boot to select sample
  fit <- lm(formula, data = d.resampled)
  return(coef(fit))
}

# regression specification
outcome <- "score"
covariates <- c("st.ratio", "lunch", "income", "english")

reg_spec <- as.formula(
  paste(outcome, 
        paste(covariates, collapse = " + "), 
        sep = " ~ "))
print(reg_spec)

score ~ st.ratio + lunch + income + english

# bootstrapping with 1000 replications
boot_results <- boot(data = CASchools, statistic = beta_hat, 
                                        formula = reg_spec, R=1000)

# standard errors - bootstrap, 1 intercept + 4 regressors
# index=1: intercept; index=2: st.ratio; ...
boot_results

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = CASchools, statistic = beta_hat, R = 1000, formula = reg_spec)


Bootstrap Statistics :
       original        bias    std. error
t1* 675.6082083 -1.083597e-01  6.03556103
t2*  -0.5603890  2.180953e-03  0.24953695
t3*  -0.3963660  1.280183e-04  0.03056716
t4*   0.6749840  2.448742e-03  0.08257985
t5*  -0.1943284  5.797645e-05  0.03397513

# the bootstrap sample mean of beta_hat_st.ratio = -0.5582080
colMeans(boot_results$t)

[1] 675.4998487  -0.5582080  -0.3962380   0.6774327  -0.1942704

# thus the bias = the bootstrap sample mean - the estimate based the original sample
# -0.5582080 - (-0.5603890) = 0.002180953 
# when constructing CI, remember to do bias-correction

# histogram of beta_hat corresponding to the st.ratio variable
# index=1: intercept; index=2: st.ratio; ...
plot(boot_results, index=2)

# get 95% confidence interval of beta_hat corresponding to st.ratio
boot.ci(boot_results, conf = c(0.90, 0.95), type = "perc", index=2) #st.ratio

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot_results, conf = c(0.9, 0.95), type = "perc", 
    index = 2)

Intervals : 
Level     Percentile     
90%   (-0.9808, -0.1819 )   
95%   (-1.0714, -0.0903 )  
Calculations and Intervals on Original Scale

トップへ戻る

偽発見率の制御 (False Discovery Rate Control)

多重性の問題。回帰推定の問題を考えます。100個の係数のうち、5つだけが結果と実際に関連しているとします。The problem of multiplicity. Imagine a regression estimation problem where \(5\) out of \(100\) coefficients have a real relationship with the result.

• この場合、真の係数セットはスパース（疎）であり、多くの\(\beta_j\)がゼロであると言えます。We would say here that the true coefficient set is sparse: most of the \(\beta_j\) are zero.

• 最良のケース では、仮説検定を実行して、これら5つの真のシグナルをすべて発見したと仮定します（すなわち、偽陰性はない）。In the best-case scenario, when you run hypothesis tests, you find all of these five true signals (i.e., there are no false negatives).

• ここで、残りの\(95\)個の係数に対して、\(\alpha=0.05\)の有意水準を用いて検定を行うとします。この場合、有用でない95個の変数の\(5\%\)について誤って有意であると結論することが期待されます。これにより、最終モデルに4.75（約5）の偽の説明変数（スプリアス・リグレッサー）が含まれることになります。Let’s assume you do find all five and test the remaining \(95\) with \(\alpha=0.05\) significance cutoff. In that case, you expect that you will erroneously conclude significance for \(5\%\) of the useless \(95\) variables, leading you to include \(4.75\approx 5\) spurious regressors in the final model.

• したがって、False Discovery Proportion (FDP)（偽発見率）は次のように定義されます：Thus the false discovery propotion (FDP) \[ FDP = \frac{\# \ \text{false positives}}{\# \ \text{tests called significant}} = \frac{5}{10} = 50\% \]

• FDP（False Discovery Proportion）は、適合したモデルの性質です。これを正確に知ることはできません。つまり、誤りを犯しているかどうかは不明です。しかし、その期待値である偽発見率 (FDR: False Discovery Rate) を制御することは可能です。The FDP is a property of the fitted model. You cannot know it: either you have made mistakes or you have not. However, you can control its expectation, the false discovery rate (FDR): \[ FDR = \mathbb{E}[FDP] =\mathbb{E}\left[\frac{\# \ \text{false positives}}{\# \ \text{tests called significant}} \right]. \]

 

次のアルゴリズムを用いることで、選択した \(q\)-カットオフ（例：一般的には \(0.1\)）に対して \(FDR \leq q\) を確実に保証することができます。You can ensure that \(FDR \leq q\) for some chosen \(q\)-cutoff (e.g., \(0.1\) is commom) in light of the following algorithm.

---

ALGORITHM 3.  Benjamini-Hochberg (BH) FDR Control

For \(N\) tests, with \(p\)-values \(\{p_1, \cdots, p_N \}\) and target FDR \(q\):

• Order your \(p\)-values from smallest to largest as \(p_{(1)}, \cdots, p_{(N)}\).

• Set the \(p\)-value cutoff as \(p^*=\max\left\{p_{(k)}: p_{(k)}\leq q\frac{k}{N}\right\}\).

The rejection region is then the set of all \(p\)-values \(\leq p^*\), and this ensure that \(FDR\leq q\). \(\quad \square\)

---

このアルゴリズムは、検定間の独立性を仮定していることに注意してください。Notice that the algorithm assumes independence between tests.

FDR control in action

# multiple regressions: 7 covariates
lm.fit.7 <- CASchools %$% lm_robust(score ~ st.ratio + calworks + lunch + computer +
                                 expenditure + income + english, se_type = "HC2")
summary(lm.fit.7)$coef

                 Estimate   Std. Error    t value      Pr(>|t|)      CI Lower
(Intercept)  6.618848e+02 1.011465e+01 65.4382631 8.746942e-220  6.420021e+02
st.ratio    -2.782033e-01 3.131453e-01 -0.8884160  3.748354e-01 -8.937650e-01
calworks    -9.174277e-02 6.369565e-02 -1.4403303  1.505334e-01 -2.169518e-01
lunch       -3.709982e-01 4.470528e-02 -8.2987567  1.518469e-15 -4.588771e-01
computer     4.100852e-04 8.207297e-04  0.4996593  6.175817e-01 -1.203255e-03
expenditure  1.743849e-03 9.054684e-04  1.9259079  5.480367e-02 -3.606548e-05
income       6.190124e-01 8.463125e-02  7.3142304  1.361586e-12  4.526495e-01
english     -2.113803e-01 3.889648e-02 -5.4344319  9.427272e-08 -2.878406e-01
                 CI Upper  DF
(Intercept) 681.767588421 412
st.ratio      0.337358497 412
calworks      0.033466221 412
lunch        -0.283119321 412
computer      0.002023425 412
expenditure   0.003523763 412
income        0.785375360 412
english      -0.134919948 412

# Benjamini-Hochberg (BH) False Discovery Rate Control, 
# Taddy (2019), p. 32, Algorithm 3 ensuring FDR \leq q, say, 0.1

# extract the p-values for the nine regression coefficients (excluding the intercept)
# and then plot p-values against their rank from smallest to largest

# slope of the BH cutoff line: q/N, where N is sample size
# decision rule: find the largest p-value that is below this line
# and call it and all smaller p-values "significant, 
# so that you expect that aroud q = 0.1 (FDR) of them are false positives


# grab p-values
#-c(1,2): drop the intercept and st.ratio 
# since st.ratio is the target variable, it should remain in the model specification
pvals <- summary(lm.fit.7)$coef[-c(1,2), "Pr(>|t|)"]  

# plot them: it looks like we have some signal here
hist(pvals, xlab="p-value", main="", col="lightblue")

# cutoff function 
BH_cutoff <- function(pvals, q=0.1){
  pvals <- sort(pvals[!is.na(pvals)])
  N <- length(pvals)
  k <- rank(pvals, ties.method="min")
  alpha <- max(pvals[ pvals<= (q*k/N) ])
  
  plot(pvals, log="xy", xlab="order", main=sprintf("FDR of %g",q),
       ylab="p-value", bty="n", col=c(8,2)[(pvals<=alpha) + 1], pch=20)
  #log="xy": both axes are to be logarithmic
  #bty="n"": no box
  lines(1:N, q*(1:N)/(N+1))
  
  return(alpha)
}

# Benjamini-Hochberg (BH) False Discovery Rate Control cutoff value
BH_cutoff(pvals) # 0.05480367

[1] 0.05480367

# at 10% FDR, we get 3 `signif'
signif <- which(pvals <= 0.05480367) 

# multiple regressions: 3 covariates
ols.fit.4 <- lm_robust(score ~., 
            data=CASchools[,c("score", "st.ratio", names(signif))], se_type = "HC2")
summary(ols.fit.4)

Call:
lm_robust(formula = score ~ ., data = CASchools[, c("score", 
    "st.ratio", names(signif))], se_type = "HC2")

Standard error type:  HC2 

Coefficients:
            Estimate Std. Error t value   Pr(>|t|) CI Lower  CI Upper  DF
(Intercept) 675.6082    6.22466 108.537 9.114e-307 663.3724 687.84400 415
st.ratio     -0.5604    0.25603  -2.189  2.917e-02  -1.0637  -0.05711 415
lunch        -0.3964    0.03033 -13.066  6.385e-33  -0.4560  -0.33674 415
income        0.6750    0.08448   7.989  1.352e-14   0.5089   0.84106 415
english      -0.1943    0.03336  -5.825  1.147e-08  -0.2599  -0.12875 415

Multiple R-squared:  0.8053 ,   Adjusted R-squared:  0.8034 
F-statistic: 463.2 on 4 and 415 DF,  p-value: < 2.2e-16

トップへ戻る

不良コントロール (Bad Controls)

回帰に職業ダミー(occupation dummy)をコントロール変数として含めるべきか？

 

以下の説明は、[Paul Hünermund氏のブログ記事]から抜粋したものです。

The following explaination is extracted from Paul Hünermund’s blogpost.

良いコントロ一ル  \(X\) は、処置変数とアウトカム変数の両方に影響を与えますが、それをコントロールすれば（それを測定することができれば）問題ありません。観察されていない交絡因子 \(U\) の存在は、この場合問題ではない。

Good Controls.  \(X\) affects both the treatment and outcome variable, but once we control for it (given that we can measure it) we’re fine – unconfoundedness holds. The presence of an unobserved confounder \(U\) doesn’t matter in this case.

 

不良コントロ一ル  \(T\)と\(X\)の間の因果関係が、今では逆の方向に向かっているということです。大学教育 \((T)\) は、将来の職業 \((X)\) に影響を与え、その結果、将来の収入 \((Y)\) に影響を与えます。もし\(X\)をコントロールすると、問題が発生します。それは、\(U\)を通る因果関係の道を広げてしまう（\(T \rightarrow X \leftarrow U \rightarrow Y\), なぜなら、\(X\)はその道のコライダーだから）、バイアスがかかってしまうのです。

Bad Controls.  Only a tiny detail has changed, namely that the causal link between \(T\) and \(X\) now goes in the other direction. College education \((T)\) affects your future occupation \((X)\), which in turn affects your future income \((Y)\). Suddenly, controlling for \(X\) causes problems. It opens up the causal path going through \(U\) (\(T \rightarrow X \leftarrow U \rightarrow Y\), because \(X\) is a collider on that path) and produces bias.

リーディング (Readings)

1. Angrist and Pischke (2015). Mastering ’Metrics: The Path from Cause to Effect. Princeton University Press.

2. 西山慶彦、新谷元嗣、川口大司、奥井亮 (2019)，『計量経済学』、有斐閣:

• 第4章 - 線形単回帰モデルの推定と検定 (pages 103 – 139)

• 第5章 - 重回帰モデルの推定と検定 (pages 141 – 207)

実証練習問題 (Empirical Exercises)

[E1]

GoogleClassroom には、デ一タファイル TeachingRatings があります。 このデ一タはテキサス大学オ一スティン校の Daniel Hamermesh 教授から提供されたもので、Amy Parker 氏との共同研究 “Beauty in the Classroom: Instructor’ Pulchritude and Putative Pedagogical Productivity,” Economics of Education Review, August 2005, 24(4): pp. 369-376 で使用されたものです。

そこには、テキサス大学オ一スティン校での 463 の講義に関する授業評価、授業の特徴、教授の特徴が含まれています。 デ一タの詳細は、TeachingRatings_Description に書かれています。 このデ一タセットの特徴の一つは、教授の「(外見上の) 美しさ」が 6 人の学生からなるパネルによって評価され、指標化されている点です。 この練習問題では、授業評価が教授の「美しさ」とどのように関係しているかを調べます。

The data file TeachingRatings contains data on course evaluations, course characteristics, and professor characteristics for 463 courses at the University of Texas at Austin. A detailed description is given in TeachingRatings_Description. One of the characteristics is an index of the professor’s “beauty” as rated by a panel of six judgers. In this exercise you will investigate how course evaluations are related to the professor’s beauty.

1. 教授の美しさ Beauty と 平均授業評価 Course_Eval との散布図を作成しなさい。 2つの変間に関係があるように見えるでしょうか。  –  Construct a scatterplot of average course evaluations Course_Eval on the professor’s beauty Beauty. Does there appear to be a relationship between variables?

2. 平均授業評価 Course_Eval を教授の美しさ Beauty で回帰しなさい。 定数項の推定値、傾きの推定値はそれぞれいくつですか。 なぜ定数項の推定値は Course_Eval の標本平均と等しいのでしょうか。 [ヒント: Beauty の標本平均はいくですか。]  –  Run a regression of average course evaluations Course_Eval on the professor’s beauty Beauty. What is the estimated intercept? What is the estimated slope? Explain why the estimated intercept is equal to the sample mean of Course_Eval.

3. 陳教授の Beauty の値は平均値であり、森教授の Beauty の値は、平均値より 1 標準偏差高い値です。 森教授と陳教授の授業評価を予測しなさい。  –  Professor Chen has an average value of Beauty, while Proessor Mori’s value of Beauty is one standard deviation above the average. Predict Professor Mori’s and Professor Chen’s course evaluations.

4. 回帰係数の大きさに関してコメントしなさい、 Beauty と Course_Eval に及ぼす推定された効果は大きいですか、小さいですか、あなたが「大きい」あるいは「小さい」と評価する基準は何ですか。 この推定された効果が因果効果なのかどうかを説明してください。  –  Comment on the size of the regression’s slope. Is the estimated effect of Beauty and Course_Eval large or small? Explain what you mean by “large” and “small”. Is the estimated effect a causal effect? Explain.

5. Beauty は授業間の評価の分散の大部分を説明していますか。  –   Does Beauty explain a large fraction of the variance in evaluations across courses? Explain.

6. 授業の種類や教授の特徴をコントロ一ルするため、変数を追加して、Course_Eval を Beauty で回帰しなさい。 具体的には、コントロ一ル変数に Intro、OneCredit、Female、Minority、NNEnglish を追加しなさい。 Beauty が Course_Eval へ及ぼす効果の推定値はいくつですか。 (2)の単回帰には、除外された変数のバイアスが含まれているでしょうか。  –  In order to control for course characteristics and professor characteristics, we regress Course_Eval on Beauty with several control variables in the regression specification including Intro、OneCredit、Female、Minority and NNEnglish. Based on this multiple regression, what is the effect of Beauty on Course_Eval? Does the simple regression in part (2) include the omitted variable bias?

7. スミス教授は平均的な容姿の黒人男性で、英語を母国語としています。 彼は 3 単位もの上級科目を教えています。 スミス教授の授業評価を予測しなさい。  –   Professor Smith is an average looking black male whose native language is English. He teaches 3 credits of advanced courses. Predict Professor Smith’s class evaluations.

[E2]

ここでは、アメリカにおける男女間の賃金格差を題材としたケーススタディを検討します。 Codebook.txt には労働者の職業関連特性の説明が記載されています。 Pay.discrimination.Rdataは、賃金および経験（exp）、性別、教育などの職業関連特性に関するCPS（2012年）のデータです。Here we consider the case study of the gender wage gap in the United States. Codebook.txt contains the description of worker job-relevant characteristics. Pay.discrimination.Rdata: the CPS (2012) data on wages and job-relevant worker characteristics, such as experience (exp), gender, education.

1. 同じ特性を持つ男性と女性の間で、賃金にどのような違いがあるでしょうか？回帰モデルの仕様をどのように決定したか説明してください。What is the difference in wages between men and women with the same characteristics? Explain how you determined the regression specification.

2. 回帰推定量の異質分散性に頑健な標準誤差を推定してください。また、各回帰係数に対する仮説検定を実施してください。Estimate the heteroskedasticity-robust standard error of your regression estimator. Conduct hypothesis testing on each regression coefficient.

3. ブートストラップを使用して、回帰推定量の標準誤差を推定してください。Estimate the standard error of your regression estimator by using bootstrapping.

トップへ戻る